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Abstract. The evolution of rotating, isolated clusters of stars up to 
core-collapse is investigated with n-body numerical codes. The simula- 
'-^ | tions start off from axisymmetric generalisations of King profiles, with 

|"y ^ . added global angular momentum. In this contribution we report on re- 

sults obtained for two sets of single-mass cluster simulations. These con- 
firm the more rapid evolution of even mildly-rotating clusters. A model 
is presented with rotational energy comparable to u>Centauri's; it reaches 
core-collapse in less than half the time required for non-rotating clusters. 
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1. Background 

o 

Star clusters are self-gravitating Newtonian systems of choice where to brew 
complex gravitational dynamics (Meylan & Heggie 1997 for a review). Observa- 
tions of old, globular, stellar clusters have led to the formulation of spherically 
symmetric dynamical models of equilibria. The most successful and universally 
studied one-integral spherical models are the King (1966) profiles. However, 
surveys of up to 100 Milky Way clusters have found small but significant depar- 
tures from spherical symmetry (White & Shawl 1987): fits to their projected 
isophotes yield ellipticities <e>= 1— <a/b>m 0.07 ± 0.01. A study of 173 
clusters in M31 found <e>= 0.09 ± 0.04 (Staneva, Spassova & Golev 1996). 
Observations of young clusters in the Large Magellanic Cloud revealed isopho- 
tal contours with ellipticities as large as e = 0.3 (Elson, Fall & Freeman 1987; 
Kontizas et al. 1990). This raises the possibility that clusters are formed as 
strongly flattened structures which then evolve towards rounder configurations 
(cf. Frenk & Fall 1982; Boily, Clarke & Murray 1999; Theis & Spurzem 1999), 
and brings up important theoretical issues concerning processes which may drive 
this evolution. 

Rotation stretches any stellar association along a preferred axis: observa- 
tions of the clusters wCentauri and M13 have shown that they are flattened by 
rotation (Meylan & Mayor 1986; Merritt, Meylan & Mayor 1997; Lupton, Gunn 
& Griffin 1987). Thus angular momentum, measured or possibly lost during 
evolution, offers a way to account for the morphology of clusters. Yet to date 
there are few evolutionary models of clusters with initial angular momentum. 
We have started on a project to develop three-dimensional dynamical models of 
rotating star clusters. In this articles results for two n-body models of isolated 
clusters are presented. Previous modelling of rotating clusters is reviewed first. 
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2. Gas and Fokker-Planck Models of rotating clusters 

Agekian (1958) considered the effects of angular momentum diffusion on the 
equilibria of rotating fluid masses of uniform density. In his analysis, concentric 
spheroids rotating about their minor axis become rounder in time when the 
spheroids have initially an ellipticity e = 1 — a/b < 0.735, where a and b are 
the minor and major axes. Shapiro & Marchant (1976) integrated the equations 
of motion for this fluid in the limit of adiabatic (slow) diffusion of momentum. 
Angular momentum losses are driven by mass elements moving in the direction of 
the stream leaving the system at a rate higher than those moving in the opposite 
direction. The energy required for escape comes from 'heat', attributed to two- 
body encounters. Thus, angular momentum losses are accrued over a local 
two-body relaxation timescale, t co \, which is inversely proportional to the mass 
density (oc In practice this hinders applications of the results to actual 

clusters, which show centrally peaked density profiles (Meylan & Heggie 1997). 
Nevertheless, the framework set by Agekian provides a start in linking rotating 
bodies and observed (non-rotating) globular clusters. 

With zero rotation, the central region of a cluster evolves towards a cusp 
in density during what is known as the gravothermal catastrophe. Does rota- 
tion stop the formation of a cusp? Hachisu (1979, 1982) discussed the time- 
evolution of self-gravitating cylindrical distributions of gas with angular mo- 
mentum. He predicted a runaway collapse of the central region whenever an- 
gular momentum is expelled faster than a critical rate (see also Lagoute &; 
Longaretti 1996). Hachisu dubbed this the 'gravo-gyro catastrophe', by anal- 
ogy with the non-rotating case. These were until recently the only evolution- 
ary models of rotating clusters. Two-dimensional orbit-averaged Fokker-Planck 
methods have now also been developed to address this issue. Following Good- 
man's (1983) approach, Einsel Sz Spurzem (1999) integrated the Fokker-Planck 
equation in energy-momentum space [E, J z \. Their initial configurations are 
truncated King models with added bulk motion. This velocity field takes the 
form of a Maxwellian distribution such that the mean velocity scales in propor- 
tion to radius away from the centre, then drops off at large radii. Their adopted 
axisymmetric distribution function (cf. Lupton, Gunn & Griffin 1987) 



where (3 is the inverse square central velocity dispersion and f2„ an angular veloc- 
ity. The initial conditions are fixed by specifying the dimensionless parameters 



ie, the scales of angular momentum and gravitational potential, respectively. 
The latter is the King parameter. 

In the 2D Fokker-Planck models core-collapse proceeds on a much shorter 
timescale than in the non-rotating case, confirming Hachisu's early intuition. 
However the central angular velocity does not increase at the high rates expected 
during the on-set of a gravo-gyro catastrophe; however near the end of core- 
collapse the central velocity dispersion bears the same relation to the central 
density as in the non-rotating self-similar collapse. This leaves open the question 



f(E, J z ) oc exp (-pn o J z ) ■ [exp (-0E) - 1] 



(1) 
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of what controls the final phase of evolution in these systems, ie whether or not 
rotation truly survives up to core-collapse. We chose to approach this problem 
using three-dimensional numerical integration; the setup is summarised below, 
followed by results and a discussion. 



3. Basic properties 

Self-consistent n-body realisations of the distribution function were obtained 
from the equilibrium Fokker-Planck code FOPAX developed by Christian Einsel. 
The models are fully specified once values are assigned to (W ,uj ). Figure [j] 
illustrates the properties of a set of models of 10,000 particles with W = 6.0 and 
four values of oj . The model clusters rotate about the z-axis and the equator lies 
in the x-y plane of a Cartesian coordinate system. Rotation causes the cluster 
in equilibrium to flatten down the z-axis and this is shown from computing the 
components of the inertia tensor J^- for a series of twenty concentric spherical 
shells of equal mass dM. We define 



rj[rk] = 1 — — , V particles in — dr < r < + dr . (2) 

Ixx ~\~ lyy 

The parameter rj = when the mass within a shell is distributed isotropically; 
7] < (or, > 0) when the distribution is anisotropic oblate (or, prolate). For the 
spherical model uj q = we found indeed near-zero values of rj at all radii. 

Models with rotation have lo q ^ and a range of values for r\ increasing with 
it. Note that all models, save one with u) Q = 0.8, have values of rj compatible 
with sphericity at the centre. At larger radii, the models are all distinguished 
from one another. 

Fast-rotating models need be more compact in order to sustain the accrued 
centrifugal force, 

v 2 

centrifugal force = — = rVL 2 , 
r 

which must always be smaller than the gravitational force, giving the condition 

n> < °£ . (3) 

' s 

Thus at constant mass M the system radius r s must be smaller to allow for 
larger angular speed f2. This is illustrated on figure l], which displays r] and f2 
computed from the same set of particles. All models show fi decreasing with 
radius. Note that the curves are consistent with solid-body rotation in the 
core-region. Further out f2 declines to near-zero, in a trend opposite that of rj. 
The core remains roundish despite the large angular speed because the gravity 
is relatively stronger there than near the edge, and so random motion of the 
particles dominate over streaming motion. Overall the fraction of kinetic energy 
invested in streaming motion ranges from 0% to 4%, 14% and 26% in increasing 
order of u . For comparisons, the cluster wCentauri invests perhaps as much as 
22 % of its kinetic energy in rotation (Merritt, Meylan & Mayor 1997). 
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See GIF figures attached 



Figure 1. Initial profile of models with W Q = 6.0 and four dif- 
ferent values of u . The parameter r\ defined in (|2|) is a measure of 
anisotropy, while f2 is the angular speed at radius r (in model units). 
Both quantities are averaged over spherical shells. 

4. N-body simulations 

The code NBODY6++ is an Aarseth-type integration code based on a Hermite 
expansion of the variables in time (Aarseth 1999). It has been ported to parallel 
architecture (Spurzem 2000); the calculations were performed on CRAY com- 
puters linked up with MPI library. The code treats particles as point-masses and 
stellar evolution options were switched off. The chain-regularisation algorithm 
for hierarchical stellar encounters as well as the standard 'KS' regularisation 
(Mikkola & Aarseth 1998; Aarseth 1999) ensures high-precision integration dur- 
ing close interactions. Only simulations with N = 5,000 equal-mass particles 
will be discussed. There are no external tides. 

Figures || & [3| illustrate the time-evolution of the models. The central 
density, total angular momentum and mean and core radii are plotted as function 
of time in units of the two-body relaxation time t co \ (see Meylan & Heggie 
1997; Casertano &: Hut 1985). (Note: the lengths were normalised to their 
initial values.) The top panels show evolution for the case of uj = 0.5, the 
bottom set for co = 0.8. Looking at these diagrams we find an evolution of the 
central density similar to the standard case with no rotation: the contraction 
of the central region leads to more close encounters and ejection, hence further 
contraction ensues, etc, until t ~ 5 t co \ when the density peaks sharply, indicating 
core-collapse: at the end of the simulations r c as 0.03 and 0.015, respectively, for 
<jj = 0.5 and 0.8. At constant energy, core-contraction drives the expansion of 
the outer envelope and hence the mean radius expands rapidly at core-collapse. 
Note a subtle but noticeable difference between the two simulations, namely 
that the cluster with initially more rotation evolves faster; this is particularly 
visible in a comparison of radii at fixed time. A more convincing demonstration 
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Figure 2. Time-evolution of the central density (left-most pan- 
els) , total z-angular momentum (middle) and the mean- and core-radii 
(right-hand panels) for runs with W Q = 6.0. 

of the fast evolution of such clusters follows if we recall that in this unit of time, 
clusters without rotation reach core-collapse in around 12 t co \, which is more 
than twice as long. 

To appreciate how many stars might be lost to galactic tides, were a tidal 
field present, we imagine the cluster orbiting the galaxy on a circular orbit. A 
tidal radius may then be defined from the initial configuration, by computing the 
radius rt at which the mean density at time t equates the initial mean density: 
nit) =<r> [0] x (M^J/MfO]) 1 / 3 . If we label as escapers all stars found outside 
2 rt(t), we obtain an estimate of the number of stars likely to leave the cluster 
on a timescale short compared with t co \; the angular momentum they carry with 
them is deduced from summing up all the momenta of the stars left behind, 
and comparing with the initial value. Implementing this algorithm, we found 
the run with uj = 0.5 (top panels) would have lost 3.7% of the initial angular 
momentum over the time of evolution, but only 1.0% (51:5000) of its mass. The 
second model, with more rotation, would have lost 0.96% (48:5000) of its mass, 
but only 1.4% of its total momentum to such escapers. This shows how the 
cluster redistributes angular momentum efficiently within itself, such that the 
core evolves towards core-collapse despite added rotational support. Evolution 
in the core is faster, the faster the core rotates initially, since the cluster is more 
compact (cf. Eq. ||| and Fig. [j]) which speeds up two-body effects. 

Figure || graphs the specific angular momentum of the uj q = 0.8 cluster as a 
function of radius r for three different times. For comparison, two components 
are given: the z-axis component about which the cluster rotates; and the x- 
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Figure 3. Distribution of specific angular momentum components 
for three time-slots. See text for explanation of the black squares. 

axis component. Dividing the cluster in ten concentric shells, we computed 
L = r x v and summed up the momenta in each shell: the result is the series 
of black squares shown on the figure. Initially the (net) z-angular momentum 
increases steadily from the centre, outwards; the symmetry of the figure would 
make the sum over L^'s cancel out, and the black squares have been left out for 
this quantity. Evolution is monotonic, with the net z-momenta inside r = 2pc 
decreasing, from which we deduce that an increasing fraction of the momentum 
is transferred to the volume > 2pc. Notice on figure || that the stars form a 
core around r = 0.5pc in the final stage of the simulation (right-most panels). 
It is not clear whether this is the result of an m = 1 (lopsided) instability, 
attributable to the dynamics of the system (from a d.f. point of view), or a case 
of core-wandering, likely due to the small number of particles inside r = 0.5pc 
(Sweatman 1993). 



5. Conclusion 

The faster evolution of clusters with rotation has been illustrated with two sam- 
ple runs. The time to core-collapse we found from three-dimensional n-body 
simulations are in agreement with two-dimensional Fokker-Planck calculations 
(Einsel & Spurzem 1999): the collapse time of 5.4 t co \ obtained for the u> = 0.8 
agrees with the Fokker-Planck solution of 5.6 t co \ for these parameters. This in- 
creases confidence in the results up to core collapse, obtained with two different 
algorithms. 
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Discussion 

A.Eckart: What happens to the rotation of the cluster after core-collapse? 

CM. Boily: I stopped my simulations precisely at core-collapse. What you 
measure in the envelope depends a lot on e.g. the galactic tide, which is absent 
in these simulations. In the core region I expect evolution will proceed much 
as in the standard case without rotation, since the cusp is isotropic and carries 
little net momentum. More detailed modelling is needed here. 
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